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I. INTRODUCTION 
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We analyze non-central heavy- ion collisions at the relativistic energy within a full (3-f 1) dimen- 
sional hydrodynamic model. First, the initial parameters in the hydrodynamic model are chosen so 
that we reproduce the experimental data of both the rapidity distribution and the slope of transverse 
mass distribution in central Pb 4- Pb collisions at 158 A GeV at the CERN SPS. We next study 
the validity of the ansatz of wounded nucleon scaling for the initial condition of the hydrodynamic 
model through analysis of the rapidity distribution of negative pions in non-central collisions. More- 
■ over, a simple Monte Carlo method is introduced to exactly calculate the particle distribution from 

I resonance decays. 
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04 , The main goals in the physics of relativistic heavy-ion collisions are not only the discovery of a new state of de- 
confined nuclear matter, the quark-gluon plasma (QGP), but also the investigation of thermodynamical aspects of its 
new phase, i.e., the equation of state (EOS), the order of phase transition between the QGP phase and the hadron 
phase, or the critical temperature [Q. The hydrodynamic model is appropriate for this purpose since thermody- 
namical variables such as temperature, energy density, and pressure are basic ingredients of the model. Hydrody- 
, namics is directly connected with the EOS of nuclear matter, thereby we describe the phase transition phenomena 
in a straightforward manner in terms of hydrodynamics; on the other hand, it is somewhat difficult to discuss the 
confinement-deconfinement phase transition by microscopic kinetic theories. Therefore, if hydrodynamics is applicable 
' to the description of space-time evolution of nuclear matter, it has some advantages over microscopic kinetic theories 
[ when we investigate the thermodynamics of the QGP. 

In addition, the hydrodynamic model enables us to analyze the collective flow of nuclear matter produced in heavy- 
ion collisions. In the central Pb + Pb collisions at 158 A GcV at the CERN Super Proton Synchrotron (SPS), 
the inverse slope parameters of the transverse mass spectra for the non-multistrange hadrons tt, K, p, and d are 
parametrized by the two common values, the freeze-out temperature and the radial velocity This implies that 
(-H , these particles constitute the radial flow, and that the local thermalization among those particles is achieved in 
central collisions. The experimental results suggest the hydrodynamic description for the expansion of nuclear matter 
seems to be suitable in central collisions at the SPS energy. On the other hand, it is an open question whether the 
hydrodynamic picture is valid even in non-central collisions. 
^ ' As well as central heavy-ion collisions, non-central events contain many attractive phenomena; these are important 
in understanding the phase transition between the QGP phase and the hadron phase. From the hydrodynamic 
point of view, there are two interesting topics in non-central collisions. One is the elliptic flow, which is one of the 
anisotropic collective flows on the transverse plane. Fourier analysis of azimuthal distribution for measured particles 
is often used to discuss the anisotropic transverse flow quantitatively [|| . The first and the second Fourier coefficients 
stand for the directed and the elliptic fiow respectively. OUitrault showed that the elliptic anisotropy is largely 
affected by the EOS of nuclear matter [Q. Sorge also showed that elliptic fiow is highly sensitive to pressure at the 
early stage of heavy- ion collisions Q . Heiselberg and Levy discussed that the centrality dependence of elliptic fiow 
may have an unusual structure when the matter at the early stage crosses the phase transition region as varying 
the impact parameter joj. Voloshin and Poskanzer showed that the centrality dependence of elliptic fiow is a good 
indicator for the thermalization of nuclear matter 0. Experimentally, the elliptic flow was observed at CERN SPS 
1^,^. The NA49 Collaboration found that the second Fourier coefficient of azimuthal distribution for charged pions 
is about 4 % around the midrapidity and the low transverse momentum [||. The WA98 Collaboration also obtained 
a similar result for positive pions. In addition to this result, they observed a remarkable phenomenon for positive 
kaons; only positive kaons are emitted out of the reaction plane at the SPS energy B. The other topic in non- 
central collisions is the "nutcrack" phenomenon advocated by Teaney and Shuryak [[lO|P According to their results, 
the unusual structure of nuclear matter can appear if the QCD phase transition occurs in non-central collisions 
at Relativistic Heavy-Ion Collider (RHIC) energies. Along the lines of these works, Kolb et al. studied these two 
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topics by using hydrodynamic simulations |Tl|. The theoretical studies mentioned above are, however, all based 
on the (2+1) dimensional hydrodynamic model P, p^|Jll|| or on microscopic kinetic theories Therefore, full 

(3+1) dimensional hydrodynamic simulations are indispensable for a comprehensive understanding of the expansion 
stage of heavy-ion collisions and the phase transition of nuclear matter in the full phase space. Here we study non- 
central heavy-ion collisions by using a fully three-dimensional rclativistic hydrodynamic model. In this paper, we 
compare the numerical results of the rapidity distribution with the experimental data in non-central collisions at 
the SPS energy, before discussing each topic mentioned above. This work is the first analysis by means of a full 
(3+1) dimensional hydrodynamic model. From a hydrodynamic standpoint, it is very important to choose the initial 
condition to reproduce the single particle spectra of hadrons since the predictive power of hydrodynamics is rather 
limited. The hydrodynamic analyses of flow, dilepton spectra, HBT, and so on, are really meaningful only after tuning 
the hydrodynamic model. 

The paper is organized as follows: In Sec. ||, we briefly review the relativistic hydrodynamic model. The rclativistic 
hydrodynamic equations with a relevant EOS are introduced in this section. We also discuss the initial condition 
of the hydrodynamic model to describe the non-central heavy-ion collisions. In Sec. Ill, we show how to calculate 
the particle distribution within the hydrodynamic model. In particular, we focus on the contribution from resonance 
decays to the particle spectra. In Sec. pv| , we show the numerical results of hadron spectra in central and non-central 
collisions. We examine in this section whether the ansatz of wounded nucleon scaling, which is widely used in the 
two-dimensional hydrodynamic model [^ UlO|JTl[ |, is really valid for the initial condition. Our results are summarized 
in Sec. 0. In the Appendix, we briefly show how to obtain the contribution from resonance decays by using a Monte 
Carlo method. 



II. RELATIVISTIC HYDRODYNAMIC MODEL 



A. Relativistic hydrodynamic equation 



Assuming local thermal equilibrium for nuclear matter produced in heavy-ion collisions, we describe its space-time 
evolution by a relativistic hydrodynamic model. Relativistic hydrodynamic equations for a perfect fluid represent 
energy and momentum conservations 



d^T^'^ix) = 0, 
T'"'{x) = [E{x) 



P{x)]u^'{x)u''{x)- P{x)g^"', 



and baryon density conservation 



dp.n'^ix) = 0, 

n'^{x) = nB(x)u^(x), 



(1) 
(2) 

(3) 
(4) 

where E, P, ub, and = (7,7V) are, respectively, energy density, pressure, baryon density, and local four velocity. 

When one analyzes central heavy-ion collisions, one can impose the cylindrical symmetry along the collision axis 
on the hydrodynamic equations. The numerical simulations of the hydrodynamic model, with cylindrical symmetry, 
are relatively easy to perform due to the reduction of spatial dimension. Hence some groups |12|-p^ have studied the 
hadron spectra in central collisions at the SPS energies by their own hydrodynamic models without using Bjorken's 
scaling solution |]l7| which helps to simplify the solution. On the other hand, in order to analyze non-central collisions, 
and moreover, to obtain the rapidity dependence of physical quantities, we numerically solve the hydrodynamic 
equations without assuming cylindrical symmetry along the collision axis or Bjorken's scaling solution [ [L7[ . We adopt 
Cartesian coordinates and rewrite the hydrodynamic equations (|l|) and (|^) as follows: 



dt 



/C/i 
U2 
U3 

\ U5 



dyP * 

d,P 
VPv 




= 0, 



(5) 



where 



/C/i\ 
U2 ^ 

U4 



/ -f^{E + P)v, \ 

' J^{E + P)Vy ^ 

-l\E + P)v, 
-f^iE + P)-P 
\ iriB / 



(6) 
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As shown in Fig. we can define the direction of x-axis as the impact parameter vector and that of z-axis as the 
colhsion axis. There are two spatial symmetries with respect to the transformations y —^ —y and/or {x, z) — > (— —z) 
for the system even in non-central collisions]^ Because of the above symmetries, we numerically solve Eq. only in 
the region — a;niax < < a^max, < y < j/max, and Q < z < Zmax- We take a;max and ?/max to be sufficiently larger than 
the nuclear diameter. The longitudinal size in a numerical simulation depends on the lifetime of the fluid: Zmax > tf- 
These spatial symmetries considerably reduce the computational time of hydrodynamic simulations. If needed, it is 
easy to obtain the solutions in the other region y < Q and/or z < by the above transformations. 

The algorithm we used to solve the hydrodynamic equations is the Piecewise Parabolic Method (PPM) This 
method is a kind of Godunov scheme. The PPM is known as a very robust scheme for solving a non-relativistic 
gas equation with a shock wave. We have extended the PPM scheme of Eulerian hydrodynamics to relativistic 
hydrodynamic equations. This enables us to describe precisely the space-time evolution of nuclear matter with phase 
transition. It should be noted that the PPM is a higher order extension of the piecewise linear method such as the 
rHLLE algorithm 

B. Equation of state 

In order to solve the hydrodynamic equations, we need to specify the equation of state, especially, pressure as a 
function of energy density and baryon density P{E,nB). Since the hydrodynamic flow is accelerated by the pressure 
gradient, the space-time evolution of fluid is governed by the model EOS. If the QGP phase appears in nuclear 
collisions, the quark matter expands, cools down, and crosses the softest point [Q. The acceleration of a fluid 
element is considerably suppressed near the softest point since the ratio of the pressure to the energy density takes 
its minimum as a function of the energy density at this point. Consequently, the hydrodynamic flow at freeze-out is 
expected to be reduced. The softest point usually lies on the boundary between the QGP phase and the hadron phase 
on the T-jj. plane. Therefore it is expected that the hadronic observables concerning the hydrodynamic flow contain 
information about the EOS of nuclear matter near the phase transition region fl^, po|j2l| ] . 

We numerically solve Eq. (^) under a relevant EOS and obtain E, ng, and Vi at each space-time point from the 
numerical solutions Ui by using P{E,nB) and Eq. (|). For an EOS, we adopt the same model as the one in Ref. [E2|. 
The model shows strong first order phase transition between the QGP phase and the hadron phase. We assume, for 
simplicity, the QGP phase consists of massless and free u, d, s quarks, and gluons. On the other hand, we adopt the 
resonance gas model for the hadron phase. The hadron phase is assumed to be composed of mesons and non-strange 
baryons up to the mass 2 GeV |p3| . It is well known that pressure in the resonance gas model overcomes that of 
the QGP phase in high baryon density and low temperature region p^ . To avoid this problem, we also adopt an 
excluded volume correction to suppress the pressure in the high baryon density region. We construct a first 
order phase transition model by the usual technique of Gibbs' phase equilibrium condition, Tqgp — Thadron (thermal 
equilibrium), /iqcp — A'hadron (chemical equilibrium), and Pqgp = -Phadron (kinetic equilibrium). A bag constant 
B^^^ = 233 MeV is chosen so that the phase transition occurs at the critical temperature 160 MeV at zero baryon 
density. The resultant pressure as a function of energy density and baryon density is represented in Fig. ||. Due to 
the Fermi energy of baryons, there exists a minimum energy density E^i^inB) at a fixed baryon density. For an 
illustrative purpose, we redefine the pressure as zero in the low energy density and high baryon density region in 
Fig. 1^. The latent heat AE(nB = 0) in our first order phase transition model is 1.4 GeV/fm'^. 

The nuclear matter produced in Pb + Pb collisions may be isospin asymmetric since the lead nucleus is composed 
of 126 neutrons and 82 protons. As well as other works based on the hydrodynamic model, we assume that our model 
EOS can be used for the isospin asymmetric matter produced in Pb -|- Pb collisions. 

C. Initial conditions 

It is difficult to describe initial hard NN collisions and pre-thermalization of the excited matter by using a one- 
fluid model, since hydrodynamics for a perfect fluid is applicable only to equilibrated mattern Hence we apply the 



'^The latter symmetry exists only in symmetric collisions such as Pb + Pb or Au + Au. 

■^It should be noted that Brachmann et al. studied, within the f/iree-fluids model, the whole of the stage of heavy-ion collisions 
including the compression of nuclear matter [Q. By using this model, they discussed entropy production at the compression 
stage. 
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hydrodynamic model after the system is supposed to be thermalized. We can choose the initial time to when the 
incident Lorentz contracted nucleus finishes penetrating the target nucleus in the center of mass system. In the 
Pb+Pb 158 A GeV collision, we obtain to = zq/vq — 2R/{voj) ^ 1.44 fm. Here vq is the velocity of incident nucleus 
in the center of mass system, i? is a radius of lead nucleus, and zq — 2R/j is a Lorentz contracted nuclear diameter. 

At the initial time to, we suppose the initial energy density and baryon density are factorized into the longitudinal 
and transverse parts as follows: 



E{to,x,y,z;b) 
nB{to,x,y,z;b) 



E{to,z)W{x,y;b), 
nB{tQ,z)W{x,y;b). 



(7) 
(8) 



Here 6 is a impact parameter and is a transverse profile function. We can parametrize the z dependence of the 
initial energy density as Bjorken's solution p7| : 



E{to,z) = EQ 



V4 



to 



[zq- 



(9) 



where we assume almost all the matter is in the QGP phase at to and its phase obeys the ideal gas EOS, P{E, ub) = 
i?/3. Since Bjorken's solution diverges on the light cone by definition, we introduce the additional factor 6{zo— \z\) so 
that the system has finite energy. An adjustable parameter zo{< zo) is to be chosen so that we reproduce the rapidity 
distribution of negative hadrons. We emphasize that we have employed Bjorken's solution as an initial condition. This 
is in contrast to Refs. |4|jl^Jll[], in which Bjorken's boost invariant solution was always assumed to be valid during 
expansion of nuclear matter. In Bjorken's picture, the resultant physical variables do not depend on the rapidity Y 
anymore fl^ . We also suppose that the baryon density as a function of z is represented by 



nB{to, z) = UBo 



to 



[Zo- I z 



(10) 



For the scaling solution, one takes a = 1; while we take a to be larger than unity in order to obtain the two peaks 
in the rapidity distribution of net proton, which are seen in the experimental data at the SPS energy [^7|] . Initial 
longitudinal flow is parametrized by 

(11) 



Vzito,z) = tanh j 9{zo- \ z |). 

This parametrization for the longitudinal velocity is different from Bjorken's solution: v-q-j{z) ~ z/to- As for the 
initial longitudinal velocity, this parametrization, rather than Bjorken's solution, enables us to reproduce the rapidity 
distribution well. We also assume transverse flow vanishes at the initial time: Vx{to) = i^yito) = 0. This means that 
transverse flow is produced and accelerated only by the pressure gradient perpendicular to the collision axis. 

The transverse density profile at to is supposed to be proportional to the number of wounded nucleons Nw in a 
colhsion i,|lO|,|ril: 



dNw 



(r;b) 




fTi„n(r- I) 



A 



n 




A 



where T, is the thickness function for nucleus i 



r,(r) - J dzp(^Vr^+z^) , 



[12) 



(13) 



and r is a vector on the transverse plane (see also Fig. |i|). For nuclear density distribution, we use a Woods-Saxon 
distribution 



p{r) 



Po 



exp[{r - R)/S] + r 



(14) 



where 6 — 0.5 fm, R = 1.12 x A^^"^ fm, and po = 0.159 fm ^. crin(= 30 mb) represents the NN inelastic cross section 
at y^SNN ^ 20 GeV. We define a transverse profile function W{x, y; b) which is normalized to have unity at the origin 
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in the case 6 = fm, i.e.^{x,y;b) = ^^(r; b)/^i^(0; 0). By shifting b in W{x,y;b) and without changing other 
parameters in Eqs. (|9|)-(pl|), we easily obtain the initial conditions for non-central collisions specified by the impact 
parameter b. It should be noted that the spatial distribution of initial collisions on the transverse plane is discussed 
in Rcf. and the procedure to obtain the density of interactions by shifting b is called wounded nucleon scaling. 
In this paper, we use the same form for the parametrization of transverse profile in the initial condition. The initial 
baryon density at origin n bo is chosen so that the number of wounded nucleons equals the number of baryons in the 
fluid: 

J (fr ^^^^ (r; ^) ^ J dxdydz nsito, x, y, z; 0). (15) 

Taking into account the spectators is beyond the scope of the present paper. Hence we assume that the spectators 
leave from the excited matter instantaneously and that its effects on the space-time evolution of hot matter or on the 
particle distribution are neglected. 



III. PARTICLE DISTRIBUTIONS 



A. Particle directly from freeze-out hyper-surface 



We assume that the hydrodynamic picture is valid above freeze-out energy density Ef and that the particles do not 
interact with each other below Thus we can define the freeze-out hyper-surface S and its elements da^. The 
numerical results of hydrodynamic simulation, i.e., temperature T{x), chemical potential fJ.{x), the local four velocity 
u^{x), and the volume element da^{x) on the freeze-out hyper-surface S, give us the momentum distribution through 
the Cooper-Frye formula^ [|2^ 

j^dN^^f p^^da.jx) 

d^p (27rr J^exp{[p'^u,{x)-fi{x)]/T{x)}Tl ^ ' 

Here d is the degree of freedom for the particles considered in the spectrum. By using this formula, we obtain the 
rapidity distribution 

dN f , , dN 

and the transverse mass distribution 



1 dN f dN 

d(f>E—. (18) 



mt drritdY J d^p 
We apply this formula for the particles directly emitted from the freeze-out hyper-surface 



B. Particles feeding from resonance decays 

In addition to the particles directly emitted from the freeze-out hyper-surface, resonance decays also contribute to 
the observed spectrum. Assuming that resonance particles are emitted from freeze-out hyper-surface and that these 
particles decay in the vacuum, the multiplicity of pions through two-body decay processes is given by [|30|,pl[ 



f J{pi,(t)\^R)dpid(j) d^pR 
Nb.^.x - j Bn^^x^ I dsWnis) 

dn f p^j^d(7^{x) 



(27r)3 is exp{[p^u,(x) - ii{x)]/T{x)} T l' 



(19) 



^The decay processes of resonance particles, however, exist. See the next subsection. 

*It is well known that this formula has a problem although it has been widely used to obtain the momentum distribution. When 
one applies the formula to the space-like freeze-out hyper-surface, it counts the particles reabsorbed by the fluid. Nevertheless 
we use this formula for simplicity. 
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where pi and (j) are, respectively, the longitudinal momentum and the azimuthal angle of pions. p'^ and cIr are, 
respectively, the resonance four momentum in the reference frame and the degeneracy. For boson (fermion) resonances, 
we take a — (+) sign. Bji^-^x is the branching ratio of the decay process and Wr is the Breit-Wigner type function, 
which takes account of the finiteness of the resonance width. For Wr, we adopt the form used in Ref. |30|. The 
Jacobian of the Lorcntz transformation from the resonance rest frame to an arbitrary reference frame J(j)i, (jj^V r) 
is defined by dpid(j)* = J{pi, (f);\ R)dpid(j), where the quantities with (without) * are the ones in the resonance rest 
(reference) frame. In the resonance rest frame, the momentum of pions is given by 

We apply Eq. ( |l9[ ) to the multiplicity of kaons or protons which emerge from two-body resonance decays. We note 
that the <j) dependence of J{pi, (j);^ r) plays an important role in understanding the observed elliptic flow of pions 

When the resonance particle decays into three particles, i.e., u Stt, we write down the multiplicity of pions as 
follows M: 



{rriR- 



R^-^XiX2 ^ J Jipi,4';VR)dpid(j)BR^^XiX2 J dW 



2 niRy^W'^ — (mi + m2)^\/M^^ — (mi — TO2)^ 
2TTQ{mR,m,)W^ 



where 



d^PR dR r PRda^jx) 

Er {27rr J exp[p-j,u,{x)/T{x)]-V ^ ' 



dx 



Q{mR,mi) = / — (m_R + m^)^ - x\/ {niR - m^)^ - x\J x - (mi + m,2)^y^x - (mi - 1712)^. (22) 
J X 

(mi+m2)^ 

The decay processes under consideration in this paper are summarized in Table |. Here we can neglect the contribution 
from resonances which have a mass larger than that of A(1232). Due to the singularity in the Jacobian J{pi, (I);Vr) 
plf , it is difficult to evaluate the momentum distribution in the usual manner. So we use a simple Monte Carlo 
method to obtain the contribution from the resonance decays. The basic idea of the method is as follows: 

(1) For each fluid clement of freeze-out hyper-surface da^{x), we generate the ensemble of resonance particles in 
the fluid rest system which obeys the Bose (or Fermi) distribution with freeze-out temperature T{x) and chemical 
potential fJ-{x). 

(2) Each resonance particle, which constitutes the thermal distribution, is Lorentz-boosted by the hydrodynamic 
flow v(x). 

(3) We make each resonance particle decay under the exact treatment based on the relativistic decay kinematics. 

(4) We count the number of decay particles which enter a momentum window under consideration and obtain the 
spectrum. 

For details, see also Appendix. 



IV. HADRON SPECTRA 



Our strategy to obtain the hadron spectra in non-central collisions is as follows. First, we choose the initial 
parameters so as to reproduce the rapidity and transverse mass (momentum) distribution for the most central bin. 
Here the transverse profiles of energy density and baryon density are based on the wounded nucleon model. Next, we 
calculate the spectra in non-central collisions by shifting the impact parameter b and without changing the other initial 
parameters. Almost all the two-dimensional hydrodynamic calculations are based on the ansatz of wounded nucleon 
Tl| ]. Although the wounded nucleon scaling for the initial condition is widely employed, its validity has 
not been confirmed yet in terms of the hydrodynamic model. In this section, we check whether the ansatz is really 
valid for the initial condition of the hydrodynamic model by analyzing the spectra in non-central collisions. 
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A. Central collisions 



First, we choose the mitial condition so as to reproduce the spectra in central coUisions. We take the four initial 
parameters in the hydrodynamic simulations as Eq — 3.9 GeV/fm^, nso — 0.46 fm^^, zq — 1.38 fm, and a = 1.7. 
When a is fixed, ubo is automatically determined by Eq. (|l^). Freeze-out energy density is fixed as E{ = 60 MeV/fm^. 
Figure ^ represents the numerical results of rapidity distributions for negative hadrons and baryons in Pb + Pb 158 A 
GeV collisions. For each distribution, we accumulate the particle directly emitted from the freeze-out hyper-surface 
and the particles from decay processes represented in Table |. Experimental data for the 5 % most central events 
were obtained by the NA49 Collaboration |2^. By adjusting the initial parameters in the hydrodynamic model, we 
largely reproduce the experimental data for negative hadrons and net baryons in central collisions. The multiplicity of 
anti-protons in our estimate is about unity, so that we safely neglect the contribution of anti-proton within the model 
calculation.^ In this calculation, we fix the impact parameter b = 2.4 fm. When one analyzes the experimental data in 
central collisions within a hydrodynamic model, the impact parameter is usually chosen to be zero. On the other hand, 
the experimental group always averages over events with small but non-zero impact parameters to obtain the spectra 
in "central" collisions. This is the reason why we choose the impact parameter as 2.4 fm rather than zero. It should 
be noted that the 5 % most central events in Pb -|- Pb 158 A GeV collisions corresponds to the impact parameter 
range < 6 < 3.4 fm |Q. The multiplicity of protons within our model EOS is exactly the same as that of neutrons. 
On the other hand, the two multiplicities are not equal in Pb -I- Pb collisions since the real nuclear matter may be 



isospin asymmetric as mentioned in Sec. II B, Hence we analyze the rapidity distribution of net baryons rather than 
that of net protons. Figure § represents the transverse mass distribution for negative hadrons and protons. We found 
that we obtain a reasonable result for the slope of transverse mass distribution when we choose Ef = 60 MeV/fm'^. 
The corresponding mean thermal freeze-out temperature and chemical potential are, respectively, < Tf >~ 117 MeV 
and < fif >^ 323 MeV. The (mean) freeze-out temperature obtained here is consistent with others |2jJ32,p3,pffl . 



B. Non-central collisions 



The adjustment of initial condition for central collisions leads us to discuss the non-central collisions through the 
ansatz of wounded nucleon scaling. Let us proceed to discussion on non-central events in terms of our hydrodynamic 
model. The NA49 Collaboration reported the analysis of non-central events and discussed the centrality dependence 
of baryon stopping |34| . They divided the events into six centrality bins and estimated the impact parameter range 
for each centrality by using the Glauber model. When analysing these data, we fix an impact parameter in the 
hydrodynamic simulation for each centrality bin and the other initial parameters are the same as those in central 
collisions. The impact parameter range which is estimated by the experimental group, the fixed impact parameter we 
choose, the number of wounded nucleons N\y, and the maximum values of transverse profile function W{0,0;b) for 
each centrality are shown in Table ||. The transverse profile functions at ?/ = fm are also represented in Fig. |^. As 
the impact parameter increases, the maximum value of the transverse profile function decreases as expected. 

We first discuss the pion spectra in non-central collisions. Figure ^ represents the rapidity distribution for negative 
pions in non-central Pb -I- Pb 158 A GeV collisions. Our numerical results for each centrality are represented by lines. 
These results are in reasonable agreement with the experimental data near midrapidity. As for the pion spectra, we 
see from Fig. ^ that the ansatz of wounded nucleon scaling works well and seems to be a proper initial condition 
for non-central collisions in hydrodynamic simulation. We next discuss the baryon spectra. Figure ^ shows the 
rapidity distribution of baryons in non-central collisions. From Fig. we found that our initialization of baryon 
density leads to almost constant baryon stopping. On the other hand, the experimental data shows that the peak 
position in the rapidity distribution depends on centrality [|^, i.e., baryon stopping increases with centrality. This 
implies that, as for the initial baryon density, the simple ansatz of wounded nucleon scaling is not enough to describe 



the impact parameter dependence of baryon stopping. We also represent in Figs, p 
mass distribution for negative pions and baryons near midrapidity in non-central co 



and K3 the results of transverse 
lisions. We see the slope of the 



transverse mass distribution becomes a little steeper as the impact parameter increases. 



^The thermal and chemical equilibrium are assumed in the conventional hydrodynamic model. On the other hand, it is often 
said that the chemical freeze-out temperature is larger than the thermal freeze-out temperature in Pb + Pb collisions at the 
SPS energy. If we take into account the additional chemical potentials for all particles at the temperature between the 
chemical and the thermal freeze-out, the multiplicity of anti-protons may become larger than the one in the present paper. See 
also Ref. S. 
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V. SUMMARY AND DISCUSSION 



We analyzed the non-central relativistic heavy-ion collisions by using the (3-f 1) dimensional hydrodynamic model 
with a strong first order phase transition between the QGP phase and the hadron phase. By adjusting the initial 
parameters in the hydrodynamic model and choosing the impact parameter b — 2.4 fm, we reproduced rapidity 
distributions for negative hadrons and net baryons in central Pb -I- Pb 158 A GcV collisions. We found that we 
reproduce the slope of transverse mass distribution with the freeze-out energy density Ef = 60.0 MeV/fm^. The 
corresponding average thermal freeze-out temperature and chemical potential are, respectively, < Tf >~ 117 MeV 
and < /if >~ 323 MeV. The freeze-out temperature obtained in this paper is consistent with the analyses of the 
NA49 Collaboration p7| , ^ . We also reproduced the experimental data of rapidity distribution for negative pions in 
non-central Pb + Pb collisions with the initial condition based on wounded nucleon scaling. The ansatz of wounded 
nucleon scaling is reasonable for discussion on the initial condition of energy density at the SPS energy. On the other 
hand, we need to modify the initial condition of baryon density in order to explain the centrality dependence of baryon 
stopping. We note that the pion spectra are not sensitive to the modification of initial baryon density. We comment 
on the initial condition of the longitudinal flow. Recently, Csernai and Rohrich discussed the 'third flow' component 
of nuclear matter in non-central collisions [^7| . In the realistic evolution of nuclear matter, the steepest direction 
of pressure gradient is tilted from both x and z axes. This causes anti-flow near the midrapidity region. In our 
parametrization, the initial longitudinal flow is independent of transverse coordinates x and y for simplicity. Thereby 
we cannot discuss this phenomenon at present. If we take into account the x dependence of initial longitudinal flow, 
for example, the initial longitudinal flow increasing near (and decreasing far from) the spectators, we can describe the 
directed and third flow in our model. 

Reproduction of single particle spectra of hadrons is very important when we try to discuss other topics or make 
new predictions by using the hydrodynamic model. Therefore we strongly propose that the experimental groups at 
the SPS and the RHIC should analyze the single particle spectra in non-central collisions before they discuss the 
anisotropic transverse flow. The numerical results of hydrodynamic simulation obtained in the present paper are our 
starting point to discuss the elliptic flow and the nutcrack phenomenon at the SPS energy. One of the authors has 
already discussed the reduction of elliptic flow due to resonance decays by using the result of the present paper . 
We will discuss the centrality dependence of elliptic flow at the SPS energy and how large the reduction effect is in a 
forthcoming paper [s^ ]. 
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APPENDIX A: MONTE CARLO CALCULATION OF THE CONTRIBUTION FROM RESONANCE 

DECAYS 

Lorentz transformation for the momentum of a decay particle between the local rest system (starred) and the finite 
momentum system (non-starred) of a resonance particle R is 



ACKNOWLEDGMENTS 




(Al) 



We rewrite Eq. (Al) explicitly 



Pi 
Pt 



PRiF{Ph(l>), 



(A2) 



cos (j>* 



Pt (pi , 0) COS (j) - pm cos (f)RF{pi , 0) 



(A3) 



\/pt{Pu4>) +PRtP'^{Pu4>) - '2-Pt{Pu4>)PRt C0S((/) - <j)R)F{pi,(j)) ' 
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where 



A\ ^{Pi^'t') Pt{pi,(l))PRtCos{(j)-(j)n)+piPRi 
rriR mR[mR + hR) 

Here the independent variables which we choose for decay particles are the longitudinal momentum pi and the 
azimuthal angle (j). Thus the transverse momentum of a decay particle pt is written in terms of pi and (f>: 

Pt{Ph4') = 7 ^({E* + PlVRljR)vRt COS{(j) - (t)R) 

± ^J [E* +piVRaRy - {pf + m^)^l{l ~ vl, cos2(0 - ^r))^ . (A5) 

The Jacobian of the Lorentz transformation is defined by 

dp*dr = J{pi,cp- VR)dpid4; (A6) 

(A7) 



J{pi,^;yR) 



dpi dip 



The calculation of J is straightforward, so that we do not represent it here. The normalization of momentum space 
volume for a decay particle in the resonance rest frame is 

* dp*i dr _ ^ 



i. (A8) 

2p* 7o 2^ 

We always average the decay probability over the spin of resonances, so that the decay probability does not depend 
on pj* and (jf . Thus the normalization in the resonance reference frame is 

J{pi,(l)\YR)dpid(t) ^ ^ 



introduce a very simp. 



The Jacobian in Eq. (A7) has very narrow peaks when the resonance particle moves at a large velocity in the 
laboratory system pi] 



Due to this singularity, it is very difficult to integrate the Jacobian numerically. So we 
le Monte Carlo calculation to evaluate the momentum distribution from resonance decays. All 
input parameters in this calculation are the numerical results of hydrodynamic simulation, the temperature T, the 
chemical potential /.t, the three-dimensional fluid velocity v, and the element of freeze-out hyper-surface da^ on the 
freeze-out hyper-surface E. In the following discussion, we show how to obtain the rapidity distribution of negative 
pions, for simplicity, only from p-mesons. In this case, the branching ratio i3pO(-)_^7r-7r+(o) = 1- It is straightforward 
to extend this scheme to the cases for other resonances or the transverse mass (momentum) distribution. 
Step 1 : Evaluate the number of p'^ and which are emitted from or absorbed by the fc-th freeze-out hyper-surface 



element dcrj^: 



9R f d^pR I PRf.da'i^ 
{2TTf J Er eMPR<lTk)-l' 



^R ^ JJR_ / l^R,.^^k I 



The integrand does not contain the Jacobian, so that it is simple to carry out the numerical integration by a standard 
technique. It should be noted that N^^ is different from the net number of emitted p-mesons from the k-th fluid 
element. 

Step 2: Generate N random momentums P* (1 < j < N) for p-mesons which obey the distribution 



(All) 



Here, for simplicity, we omit the Breit-Wigner function. 

Step 3: For each N random momentum P* , generate random variables (0*,$*) whose ensemble is uniformly dis- 
tributed on the unit sphere. By using these random variables, we obtain an ensemble of p-meson with momentum 
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P* = {P*j,P*j,P*j) = (Pj'sin9*cos$*,P*sin6*sin<I>*,P*cos6*), which obeys the Bose- Einstein distribution in 
the fluid rest system. 

Step 4 '■ Boost P* with respect to the fluid velocity : 

p,^p; + v,7.(i?; + ^f^), (A12) 



where -jk = 1/^/1 - v^. 

Step 5: Generate N uniform random variables on the unit sphere {9*,(j)*) and obtain an ensemble of pion with 
momentum p* = {p*xj,Plj,Plj) = [p* sin 6** cos (j}*,p* sin 6* sin(/)*,p* cos 6**), where p* is given by Eq. (|o|). 
Step 6: Boost p* with respect to the resonance momentum Pji 



Step 7: If P'^da^k is positive, 



if P'l'da^k is negative, 



+ (A14) 



Y + and P'^duak is positive, then 



I P'^da^k I 

(A15) 

Here, N'^ {Nk) is to be proportional to the number of p-mesons which are emitted from (absorbed by) the fc-th fluid 
element. 

Step 8 : If the rapidity of a negative pion Yj which is evaluated from Pj enters in a rapidity window Y — <Yj< 
^ and Pjda^k 

AN+iY) ^ AN+iY) + (A16) 

if Yj also enters the above rapidity window but P^da^^k is negative, 

I Ptdank I 

AN-{Y) ^ AN-{Y) + (A17) 



E* 

Step 9: Repeat steps 7 and 8 for all N random variables. 

Step 10: Obtain the rapidity distribution of decay particles from the fc-th fluid element 

^(r) . ^AKiY)~^N^iY)^ 

dY ^ ' N+ + N- \ AY J ' 



It should be noted that the minus sign in the bracket of Eq. (A18) means the net number of emitted particles from 
the fc-th fluid element. For the normalization in Eq. ( A1S| ), we use the gross number + NjT since this number is 
positive definite. 

Step 11: Repeat the above steps from 1 to 10 for all fluid elements obtained in a numerical simulation of the 
hydrodynamic model. Summing over the contribution from all fluid elements on the freeze-out hyper-surface S, we 
obtain the rapidity distribution of tt^ which are from p decays: 

k 
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V 
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FIG. 1. Schematic view of a non-central heavy-ion collision at the initial time to- The left (right) figure represents the 
reaction (transverse) plane, b is an impact parameter vector. The gray areas represent the spectators which we neglect in the 
present paper. 
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FIG. 3. The rapidity distribution of negative hadrons and baryons in central collisions. We fix an impact parameter b = 2.4 
fm. The solid (dashed) line represents the rapidity distribution of negative hadrons (baryons). 




mt-mo GeV 

FIG. 4. The transverse mass distribution of negative hadrons and protons in central collisions. We fix an impact parameter 
b = 2.4 fm. We adjust the freeze-out energy density Ef to obtain reasonable agreement with the slope of each spectrum. The 
resultant E{ is 60.0 MeV/fm*. 
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FIG. 9. The transverse mass distributions of baryons in non-central collisions for each centrahty. 



TABLE I. Resonance decays which contribntc to the spt;ctrnni. Tht; ft'actions in the branching ratio conic horn the Cleb- 
scli-Gordau coctficicnts. 



Decay particle 


Decay channel 


Branching ratio 




p" — » 7r~7r"'" 


1.0 




p~ —>■ it~it'' 


1.0 






0.88 




K*" Tr-K+ 


1 X 1.0 




K*- ■K-K° 


fx 1.0 




A~ — » TT~n 


1.0 




A° TT-p 


i X 1.0 

1 X 1.0 


K- 


K*- K-Tr+ 




K*° K-n° 


1 X 1.0 
1 X 1.0 


P 


A° pTT- 




A+ ^ piT° 


fx 1.0 




A++ -^ P-K+ 


1.0 


n 


A+ ^ n7r+ 


1 X 1.0 




A° ^ n77° 


fx 1.0 




A~ — » n7r~ 


1.0 



TABLE II. The impact parameters in the hy<h'odynarnic siimilations. A^r represents the nuinl)er of wounded micleons. W 
is a transverse profile function. For details, see text. 



Binl Bin2 Bin3 Bin4 Bin5 Bin6 

range (fm) 0-3.4 3.4-5.4 5.4-7.4 7.4-9.1 9.1-10.2 10.2- 

b (fm) 2.4 4.4 6.4 8.2 9.5 10.6 

Nw 353 282 199 128 83 51 

W{0,0;b) 0.983 0.941 0.869 0.770 0.669 0.554 
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